function [ W, H ] = LSNMF( duomenys, R, niter, varargin)
X = duomenys;
[D, N] = size(X);

W = rand(D, R);
H = rand(R, N);
[w0 h0] = Popt(varargin, 'W0', W, 'H0', H);
a = duomenys;

tolfun = 1e-4;
tolx = 1e-4;

sqrteps = sqrt(eps);
nm = numel(a);
dnorm0 = 0;
for k=1:niter
   
        % Alternating least squares
        h = max(0, w0\a);
        w = max(0, a/h);
     % Get norm of difference and max change in factors
     
    d = a - w*h;
    dnorm = sqrt(sum(sum(d.^2))/nm);
    dw = max(max(abs(w-w0) / (sqrteps+max(max(abs(w0))))));
    dh = max(max(abs(h-h0) / (sqrteps+max(max(abs(h0))))));
    delta = max(dw,dh);
    
    % Check for convergence
    if k>1
        if delta <= tolx
            break;
        elseif dnorm <= tolfun*dnorm0
            break;
        elseif k==niter
            break
        end
    end
 
    w0 = w;
    h0 = h;
end

W = w;
H = h;


end
